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A method for engineering the behavior of populations of rhythmic elements is presented. The 
framework, which is based on phase models, allows a nonlinear time-delayed global feedback signal 
to be constructed which produces an interaction function corresponding to the desired behavior of 
the system. It is shown theoretically and confirmed in numerical simulations that a polynomial, 
delayed feedback is a versatile tool to tune synchronization patterns. Dynamical states consisting of 

QQ , one to four clusters were engineered to demonstrate the application of synchronization engineering 

f^*) . in an experimental electrochemical system. 
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Populations of interacting rhythmic components can produce complex behavior in biology [l], [2|], 
communications [3j, population dynamics Q], and chemistry [H, [H, 0, 1H . In biology, synchronization 
can be beneficial, such as in orchestrating the circadian rhythms in mammals, or pathological, such as 
in the occurrence of Parkinson's disease. We consider here the engineering of desirable states through 
the introduction of mild feedback, mild such that the behavior of the individual components is not 
substantially changed by the introduction of the external signal. In a previous paper [9], we have 
■ experimentally demonstrated a general methodology for engineering a target dynamical behavior in 
oscillator assemblies. The aim of the present paper is to describe the theory behind our methodology 
and to verify it by numerical and experimental studies. 
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I. INTRODUCTION 

Ensembles of self-sustained oscillators can spontaneously organize their collective dynamical behavior as a result of 
interaction amongelements. Examples can be found in biological [HQj chemical @, @, 0, Hi , and ecological systems [3], 
communications Q , as well as human activities (lfj| . Global behaviors such as synchronization are often responsible for 
^ — , the formation of certain beneficial biological functions, such as orchestrating the sleep/wake cycle (circadian rhythm) 
of mammals (Tlj . Conversely, pathological synchronization may induce serious problems, e.g., tremors in Parkinson's 
disease @ and abnormal sway in London's Millennium Bridge [T3 |. 

Interactions involving feedback among rhythmic elements are often associated with the formation of dynamical 
qq ' order. In many cases, feedback plays an essential role in sustaining dynamical stability and in suppressing complexity. 
In chemical systems, several types of complex synchronized behavior emerge by introducing global feedback, which 
' otherwise shows simple synchronization or chemical turbulence 0, [ID, [lj[ ■ It may be expected that feedback loops 
. i-h , among neuronal clusters contribute to the design of the complex dynamical functionality of the brain . In medical 
■ applications, heart pacemakers, deep brain pacemakers, and feedback control techniques have been proposed to 
eliminate pathological synchronization 0, [IB]. For applications which involve biological neurons, a mild control is 
desired to avoid side effects and to maintain the fundamental nature of neurons in the system. 

In this paper we present a comprehensive theory for designing the collective dynamics of a rhythmic population 
using external feedback and, as an application of its use, demonstrate the power of the methodology for creating 
various cluster states in both numerical simulations and experimental studies. 

This method has been shown to be extremely robust in engineering collective dynamical behavior in electrochemical 
experiments Our method utilizes phase modeling || (as opposed to a physiochemical modeling) to describe the 
dynamical behavior of rhythmic elements. The simplicity and analytical tractability of the phase model is exploited to 
design an optimal, delayed, nonlinear feedback signal for obtaining a desired collective behavior. The only properties 
required to construct the phase model are the waveform and the phase response function of the oscillator, which can be 
experimentally measured. Sec. II outlines the mathematics of the methodology. In Sec. Ill, the detailed description 
of the feedback design method is presented, which is further developed in Sec. IV for both harmonic and slightly 
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inharmonic waveforms. Sees. V and VI present examples of synchronization engineering using numerical studies via 
the Brussclator model and experimental studies using electrochemical oscillators, respectively. 



II. GENERAL METHODOLOGY 



Our methodology seeks to engineer a target collective dynamical behavior in a population of limit-cycle oscillators 
through feedback. A phase model is used to describe the collective behavior of a population of weakly coupled 
oscillators. For N oscillators with general interactions (which also admits interactions through a complex network), 
the phase model is given as 

^- = cj i + KY,H ij (<l> j -<t>il (1) 

3 

where <f>i (0 < <fii < 2tt) and u>i are the phase and the natural frequency of the oscillator i, K is a coupling strength 
(K > 0), and the 2-7T periodic function Hij(4>) is the (phase) interaction function (or, phase coupling function). The 
phase model is derived as the first order approximation of a coupled limit-cycle oscillator system, where the small 
quantity is the coupling intensity K Q . The interaction function £Ty (0) can be calculated from the properties of the 
limit-cycle oscillator i and physical interaction from the oscillator j to the oscillator i. The details of this derivation 
are presented in Sec. Mil Although, as presented in Sec. IIH) our proposed theory may deal with Eq. fl}, we mostly 
devote ourselves in the present paper to the case of global feedback. In such a case, the phase model is reduced to 

j 

Moreover, if the heterogeneity is sufficiently small compared to the feedback intensity K, we can treat the system as 
identical oscillators. In such a case, we may use the following phase model: 

= <" + f£*fo-*)- ( 3 ) 

3 

For simplicity, we outline our methodology in terms of Eq. ([3]) in this section. 

In the phase model ([3]), dynamical evolution of the system is predicted if the interaction function H((j>) and an 
initial condition are given. This comes from the fact that we may set u> = and K = 1 without loss of generality 
(using a rotating reference frame, 4> — cot — > cj> , and rescaling time, Kt — > t). The relationship between the shape 
of the interaction function and the global dynamical behavior of the system has been well studied. For example, 
the conditions which admit perfect synchrony, perfect desynchrony (the spray state) @, phase clustering [17j . and 
slow switching dynamics [l8l . [l^ . [20| are known. While it is possible for an interaction function to admit multiple 
attractors, it is preferable to have a single stable attractor (or at least a single dominant basin of attraction). If a 
coupled limit-cycle oscillator system has a phase interaction function which results in a single stable attractor, the 
system should exhibit the expected global behavior under general initial conditions. 

Two steps are required to engineer a desired target behavior: 

(i) find an interaction function H(<p) which uniquely stabilizes the desired collective behavior 

(ii) seek the physical feedback parameters which result in the interaction function found in (i) 

The difficulty of the step (i) depends on the desired target behavior. To illustrate the engineering methodology, we 
have selected a simple collective behavior which has been well characterized in terms of the phase model: perfect 
synchrony and balanced cluster states. Two different approaches can be utilized to optimize feedback parameters. 
The first approach simply requires knowledge of the precise interaction function to be targeted. However, there are 
many cases (such as phase clustering) in which there exists a large family of valid functions, each capable of producing 
the desired target behavior. By arbitrarily selecting one of these functions, the most effective means of generating 
the target behavior may be overlooked. Additionally, an arbitrary choice of an interaction function can substantially 
increase the difficulty of the feedback parameter optimization. Therefore, instead of targeting a precise interaction 
function, we place constraints on specific properties of the interaction function as required to generate the appropriate 
behavior. The optimum feedback parameters are selected such that the associated interaction function meets these 
constraints with the minimum feedback amplitude. 
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To engineer a target interaction function into a physical system, we introduce a feedback Kp(t) to some global 
parameter of the system with the following functional form: 



1 N 

p(*) = jvE^)' w 



s 

h{xi) = k n {xi(t - r n ) - ao}", (5) 

n=0 

where Xi(t) is an observable variable of the oscillator i at time t, ao is the time average of Xi, k n and r„ the gain and 
the delay of the n th order feedback respectively, and K and S the overall gain and the overall order of the feedback 
respectively. Our choice of function ([SJ is motivated from the fact that each feedback term yields different combinations 
of intensities of Fourier components and the feedback delay value r„ controls the ratio between the symmetric and 
antisymmetric Fourier components. In addition, the n th harmonic of the interaction function is efficiently enhanced 
by the n th order feedback. Thus, flexible and efficient design of the interaction function is possible. In Sec. IIII| we 
show that any target interaction function which is composed of S and lower Fourier components can be indeed 
produced by using the S overall order feedback. In particular, when the waveform Xi(t) is exactly harmonic, feedback 
parameter values {k n } and {t„} may be calculated analytically, as illustrated in Sec. IIVI For a general waveform, a 
numerical optimization is often required to determine the feedback parameters. 



III. THEORY OF FEEDBACK DESIGN 



We present a theory for designing external feedback signal yielding a desired phase interaction function. Because 
the extension to more complex situations is straightforward, it is suitable to start with the case where the oscillator 
1 is affected by the feedback signal composed as a function of the state variable of the oscillator 2. The dynamical 
equation for the oscillator 1 is then given as 

^M=F 1 [A 1 (t)]+KP(t), (6) 

where Aj = (xi,yi, . . .) T is the state variable of the oscillator i (i = 1,2), Fi(Ai) is a nonlinear function admitting 
limit-cycle oscillation, and P(t) is a feedback signal. The dynamical equation for the oscillator 2 is arbitrary provided 
that it produces nearly periodic dynamics. We define the observable variable to be x and the variable to be perturbed 
by the feedback to be y (these variables need not be mutually exclusive, which is the case in the numerical studies 
and experiments). Thus, 

P(f) = {0>( a;2 ),0,...,0} T , (7) 
where h(x) is Eq. ([5]). 

Because we are interested in mild engineering, we assume that the overall gain K is small such that the dynamical 
behavior of the system ([6]) can be approximated by the phase oscillator model, 

^-=u Jl +KY^H(d> 2 - ( f> 1 ). (8) 

3 

The interaction function H{4>2 — 4>i) is computed as 
1 f 2 * 

H(?h-<h.) = — Z{<j )1 + 6)h{4> 2 + 9)d6 (9) 
27T Jo 

where Z(<fii) and h(<f>2) are evaluated from single isolated oscillators 1 and 2, respectively. The function Z{(f>i) is 
referred to as the phase response function (or, the phase sensitivity function) of the oscillator 1, which is the gradient 
of the phase along the y-dircction on the limit-cycle orbit A± (<p), 

Z(M = ^ c ■ (10) 
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There are several ways to measure Z(<p) of a given oscillator . The function h{4>2) is obtained by first describing 
X2(t) as the function of the phase of the oscillator, a^f^)- Because 4>2 (t — r) = 4>2 (i) — w 2 r when the interaction is 
absent, we have xi($i{t — r)) = X2(4>2(t) — ^r). As a result, /i(<fe) assumes the form 

s 

M^) = E K{x{<t>2 ~ u 2 T n ) - a } n , (11) 
or, upon expanding x{4>) = J2i a;e~ 1 ^, as 

n 

J2 aie -il^ e il^r n I _ (12) 
i^O J 

Thus, the phase model is an autonomous system despite the existence of time delays in the original system dB]). Such 
an approximation is valid so long as Kt remains small (note that the dimension of K is inverse time) fl9l. [25|. 

For given Z^fa), the feedback parameters k n and r n yielding a target H((f>) are found in the following way. To 
simplify the problem, the functions are expanded into their Fourier series. H{(j>) = J2iHie~ l1 ^, Z(<j>) = J2i Zie~' a ^ 
and h((f>) = '^2, l hie~ %1 ^' (where Hi = H*_,,Zi = Z*_^ and hi = h*_,). Using these Fourier coefficients, we obtain the 
relation 

Hi = hiZ-u (13) 

where hi is the function of k n and T n . By solving a set of complex equations (|13p . the feedback parameters k n and r n 
can be determined. In theory, any interaction function composed of harmonic components Hi for < I < S can be 
constructed using a feedback signal with an overall order of S, provided that zi for I = 0, . . . , S is non- vanishing. 

It is important to point out that the flexibility of our engineering method is reduced for certain types of oscillators. 
For example, the Stuart-Landau oscillator (i.e., the normal form for the Hopf bifurcation) has Zi = for I > 2 0], 
forcing all higher harmonics (I > 2) in the interaction function to vanish regardless of the nonlinear terms in the 
feedback. A similar problem may occur in systems which contain special symmetry. For example, the Van der Pol 
oscillator has symmetry with respect to the center of the oscillation. This fact implies that Z{<j> + tt) = Z(<fi), i.e., 
Zi = for even I. Thus, the even Fourier components of the interaction function vanish. In these special cases, the 
methodology is limited to controlling only those harmonics of the H((f>) which do not correspond to the vanishing 
harmonics in Z(<j)). 

It is straightforward to extend the above arguments to a population of oscillators. Under the assumption that a 
parameter in each individual oscillator can be independently tuned online, we may consider 

lMl=F i [A i (t)}+K(0,Pi, ■■■)*, (14) 
where the function pi is fully generalized as 

ft(*)= EE ( 15 ) 

j n=0 

The corresponding phase model then reads Eq. (fTJ. The interaction function Hij{(f) is determined as a function 
of the physical feedback parameters {k^} and {r^ }; Any Hij{4>) can be designed by giving appropriate feedback 
parameters. The phase model (TT]) is very general and a large class of collective behavior can be engineered. 

A simple situation, which is the case in the numerical and experimental studies in Sees. IVland lVll is that a global 
parameter of the system is tuned by feedback. In such a case, we consider the phase model © by adopting the global 
feedback signal ([4J and ([5]). 



1 For example, see Ref.^for the analytical derivation, the software by Ermentrout, xppaut, for the numerical derivation, and Refs.l2ll.l22l.l23 
for the experimental derivation. Some of the methods are reviewed in Ref. |24| . 
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IV. USE OF HARMONIC SIGNALS 



An oscillator may appear to have a nearly harmonic waveform by the nature of oscillation or through the use of a 
lowpass filter. In the case of a perfectly harmonic waveform, feedback parameters can be explicitly calculated and the 
effect of each feedback term on the interaction function clearly understood. However, it is unrealistic to assume that 
an exact harmonic waveform can be obtained experimentally Therefore, the effect of a weakly inharmonic waveform 
is also examined. 



A. Harmonic Waveform 

We first assume that x(<f>) is an exact harmonic waveform with zero mean. Properly defining the origin of the phase, 
we may set 

= e _< * + c.c. (16) 

where c.c. indicates the complex conjugate (i.e., e 1 * in this case). 
Introducing <f) n = <j) — u)T n , we obtain 

n 

2#«) n = E c;;^™- 2 ™^" (17) 

m— 

where C™ a is a number of m combinations from a set of n elements. Since x((j) n ) n contains only / th harmonics (where 
I = n, n — 2, n — 4, . . .) the n th feedback term in Eq. (fT5)) only contributes to the I th (I = n, n — 2, n — 4, . . .) Fourier 
components of the interaction function. The feedback delay varies the ratio between the even and odd components 
of each harmonic. 

Combining the feedback terms using Eq. (|5|) yields an interaction function composed of Fourier components Hi for 
I < S. 

S n 

K<t>) -EE c^"- 2 ™)^ (is) 

n— m— 

or, for I > 

s 

hi = VfcnC^e^" (19) 

' * 2 

n—l 

where for convenience we define = if m is not an integer. Therefore 

S n 

H{4>) = E fc » E C? n z n . 2m e^ 2m - n ^ , (20) 

n— m— 

or, for I > 

s 

Hi = z_ i y^k n C n n+l e a ^, (21) 

' 2 

n—l 

Therefore, for a given target Hi (I = 0, . . . , S), the feedback parameters k n and r„ (n = 0, . . . , S) can explicitly 
obtained. Since H$ is determined solely by the 5 th term, the parameters k$ and t$ can be found by solving a 
complex equation obtained by setting / = S. The same process can be used for the Hg-i component, to solve for 
the parameters and ts-i- Since H$-2 is dependent on the S th and (S — 2) th terms, fcs-2 and rs_ 2 can be 

determined. Repeating this processes for each term, all feedback parameters can be calculated. 

As an illustration, the feedback parameters required to produce a Hansel-Manubille-Mato type [l8[ interaction 
function will be calculated. The target function H{4>) is selected to be 

H(<t>) = wa.(<f> - a) - r sin(20) = - V Q e~ 10 + ^e~ 2i * + c.c, (22) 
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where a and r > are the parameters of the function. Since the target interaction function has second order 
components, second order feedback is required. Therefore, equation (f2"0")) for S = 2 becomes 

H{<j>) = (fc + 2k 2 )z a + hz^e^e-^ + k 2 z_ 2 e 2luJT2 e- 2 ^ + ex.. (23) 

Comparing Eq. (f2"2")) and Eq. (|2"3")h we find one of the solutions to be 

T 1 T 

fc o = - | 1 , ki = -, 1 , fc 2 = 1 , (24) 

\z-2\ 2 |z_i| 2|z_ 2 | 

a - f - arg(2_i) f - arg(z_ 2 ) 

Ti = , r 2 = . (25) 

uj 2uj 

B. Slightly Inharmonic Waveform 

The effect of weak inharmonic components is considered using the waveform 

x{4>) = + ee l20 + 0(e 2 ) + c.c., (26) 
where e is a small complex number. Introducing <f> n = <j) — ujt„, for n > 1 yields 

n n — 1 

X (0„)" = ^ C n e i(n-2m)^ + £ ^ J e «(n-2m+l)0„ + e4 (n-2m-3)0„ ^ + ( £ 2^ 

m— m— 

and a; = 1. The n th term contributes to the I th (I = n + 1, n — 1, . . .) harmonic with order e. 5 th order feedback 
strongly enhances the harmonics hi of the interaction function for I < S. Therefore, S th order feedback is required to 
produce a target interaction function composed of harmonics / < S. Although the (S+ l) th Fourier component appears 
in the interaction function, it is expected to be very small [of the order of 0(eZs+i)] and can be safely neglected. 
Similarly, this result is also true in cases where the first order Fourier component of the waveform is dominant. 

When the waveform is strongly inharmonic, each feedback term enhances various harmonics, including higher order 
harmonics Hi (I > S) of the interaction function. In these situations, no analytical solution is possible, and the 
feedback parameters must be numerically optimized using Eq. @ or Eq. p3[) . In addition, usually, a high order 
feedback is required (i.e. large S) such that Zi for I > S are negligible. 

V. NUMERICAL STUDY: PHASE CLUSTERING 

Our methodology is numerically verified for a population of limi t-cy cle oscillators, using the Brusselator model, 
a simple two variable ODE system that admits a Hopf bifurcation [26J . The dynamical equations for a Brusselator 
population under global feedback are 

-Jr = ( B - l ) x ^ + A 2 x l + f(xi,yi) + jj^2 h ( x j), 

3=1 

^ = -Bx l -A 2 y l -f(x ll y l ) (28) 

where f(x, y) = ^x 2 + 2Axy + x 2 y. Here, it is assumed that the feedback signal is constructed from and applied to 
the variables Xi . Note that for convenience, the variables Xi and yi are transformed such that the fixed point is shifted 
to (x, y) = (0, 0). For a single uncoupled oscillator, the Hopf bifurcation occurs at B = B c = 1 + A 2 . The parameters 
of Eq. (|2"5]) were chosen to be A = 1.0 (so that B c = 2.0) and B = 2.3. The waveform x(4>) and the response functions 
Z{<p) along the x-direction 2 are displayed in Fig. [1] and their Fourier coefficients can be found in Table [T] 



2 The response function is calculated using xppaut. 
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FIG. 1: (Waveform x(<f>) and phase response function Z(<f>) of the Brusselator oscillator. The parameter values are A = 1.0 and 
B = 2.3, resulting in T = 6.43 and w = 0.977. 
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TABLE I: Brusselator (A — 1, B — 2.3). Fourier coefficients of the wave form x(<f>) and the phase response function Z(cj>). 



The phase of an oscillator fa is determined by measuring the times at which the orbit in (xj, j/j) state space crosses 
a particular point on its limit cycle. The phase of the oscillator fa(t) for t m —i <t<t m can be defined as 

fa{t)= 2 f-_l^\ (29) 

where t m is defined as the time at which the m th crossing occurs. Note that < fa(t) < 2ir. To quantify the amount 
of collective global order within the system, it is useful to introduce the order parameter R^: 



N 

,ik<t>j 

TV 



1 N 

-Y> 



(30) 



For a large population (i.e. large TV), Rk = for uniform phase distribution and Rk = 1 for a balanced k cluster state, 
in which the population is split into k equally populated point clusters distributed uniformly in phase (see Appendix 

ED- __ 

Phase clustering commonly appears in globally coupled oscillator systems [13, 1231 • In systems of identical coupled 
oscillators, these states always exist independent of the interaction function, such that the only outstanding issues 
to be addressed are the stability and the basin of attraction of the states. In this example, four parameter sets are 
created [one for each cluster state (n = 1,2,3,4)] with the following conditions: 

(i) the n cluster state is uniquely stable among the balanced cluster states 

(ii) the cluster state has high linear stability 

(iii) small amplitude feedback is preferable. 

For condition (i), it is convenient to use a target interaction function of the form lmH n > and Imi?j < for 
I ^ n (note that the symmetric parts KeHi are irrelevant to the stability of the balanced cluster states, so that we 
may arbitrarily set the symmetric parts). For such an interaction function, the maximum eigenvalue is given by 

Amax = — ^2iZ\ 2Z ImHni (see Appendix lAl for the details). We thus require Am2 x < 0. Satisfying condition (ii) 
requires that lmH n is large enough for high stability. To satisfy condition (iii), n th order feedback is used to generate 
the n th cluster state (i.e. S = n), since the n th cluster state requires n th order harmonics in its interaction function. 

A large number of interaction functions satisfies conditions (i) and (ii). Out of this family of interaction functions, 
the optimal feedback parameter set is selected such that it minimizes the cost function X)/=i 1^1 1 under the conditions 
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TABLE II: Brusselator population with global feedback. The target and resulting interaction functions, feedback parameters, 
and the resulting maximum eigenvalue for the n cluster state. 




5 10 15 20 5 10 15 20 

Kt Kt 

FIG. 2: (Color online) Engineering cluster states in the Brusselator model. Time traces of the order parameters R n and 
snapshots for (a)n = 1, (b)n = 2, (c) n = 3 and (d)n = 4 for N = 12 and K = 0.001. In each panel, the parameter set n 
diplayed in Table HT1 is used. In the insets, snapsh ots in the one-oscillator phase space at Kt — 20 are displayed. The dashed 
lines are the orbits of an oscillator. In panel (a), the initial conditions are shown (on the left side). Note that the initial 
condition is the same for all panels. 

shown in the left side of each column of Tabic [Til The table also displays the optimized feedback parameter sets 
(obtained numerically using Mathcmatica), the resulting ImHi (in the right side of each column), and the resulting 

maximum eigenvalue A max- 

Applying the optimized feedback parameter sets to Eq. ([25]) causes the system to approach the desired cluster 
states. The convergence of the system to the cluster states is illustrated in Fig. [21 using the appropriate order 
parameter R n . Several different random initial conditions were used for each parameter set and in each case the 
desired cluster state was obtained (not shown). 

It is worth noting that, in practice, the 3 and 4 cluster states are difficult to obtain unless high order feedback 
is used. For example, when only the linear term is used, the magnitude \Hi\ = \Ziai\ is very small for I > 3. This 
fact implies that the maximum eigenvalue of the n > 3 cluster states can not be large and negative. Therefore, the 
presence of noise or heterogeneity, if any, would destroy the n > 3 cluster states. 
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VI. EXPERIMENTAL STUDIES 



A. Experimental Setup 

The preceding theoretical work on synchronization engineering was experimentally tested using a population of 
electrochemical oscillators. These oscillators were created using an electrochemical cell which consisted of 64 Ni 
electrodes (99.99% pure) in a 3M H2SO4 solution, a Pt mesh counter electrode, and a Hg/Hg 2 S04/K 2 S04 (sat) 
reference electrode. The cell was enclosed in a jacketed glass vessel held at a constant temperature of 11°C. An 
EG&G potentiostat was used to adjust the circuit potential (V) of the cell, causing the nickel electrodes to undergo 
transpassive dissolution. The dissolution current of each electrode, Ij(t), was measured by zero resistance ammeters 
(ZRA). A resistor (R p ) was attached to each channel to induce oscillations in the electrode potential. A Labview 
based real time data acquisition computer was used to read the ZRA measurements, stream these measurements to 
the host machine, and apply the feedback signal to the potentiostat at a rate of 250 Hz. The current measurements 
were scaled: 

j/ (t) = A™-n (J (t) _ joffset) (31) 

The mean value of each channel (j° ffsot ) W as removed from the measurement, and the result was scaled by the 
amplitude of its oscillation (Aj) relative to the mean amplitude of the population (A mcan ). The host machine was 
used to continuously determine the offset and amplitude of each rhythmic element in the population. To calculate the 
feedback signal, the potential drop across the double layer, Xj (t), was determined from the scaled current measurements 

x j (t)=V(t)-I' j (t)R p (32) 

where V(t) is the applied voltage. The perturbation signal p(t) which was fed back to the potentiostat, V(t) = 
Vq + Kp(t), was calculated by taking the mean value of h(xj(t)) over every element in the population, 



1 N 



3=0 



h(x(t)) = - r «)" ( 34 ) 

where K is the overall feedback gain, N is the number of elements in the population, k„ is the n th polynomial feedback 
coefficient, r„ is the time delay of the nth polynomial feedback term, and S is the polynomial feedback order. 



B. Experimental Validation of Theory 

The synchronization engineering framework, as derived in Sec. IIVBI for weakly inharmonic waveforms, predicts 
that n th order feedback can enhance up to and including the n th harmonics of the interaction function for harmonic 
waveforms, and (n + l) th harmonics for weakly inharmonic waveforms. To test the range of the validity of this result, 
the experimental system was used to measure how the harmonics of an interaction function change with increasing 
feedback order. The operating voltage of the system was selected to be 1.110 V as this was found to be close enough 
to the Hopf bifurcation to produce a nearly harmonic waveform, but far enough to ensure that the periodic cycle was 
robust against external perturbations [Figs. [3jA) and[3]^B)]. 

A two oscillator system was used to measure the interaction function associated with the global feedback signal. 
This method of measurement was created by extending the work of Miyazaki and Kinoshita [28| to rhythmic systems 
under global feedback. By measuring the change in the period of the two elements as a function of their phase 
difference, the interaction function can be experimentally determined. 

While the 1 st order harmonic of the waveform accounts for 71% of its magnitude, this may not be sufficient to 
allow the 0(e) terms of Eq. (|2"T|) to be neglected. Therefore, it is expected that the (n + 1) order harmonics of 
the interaction function will be dominant. Figure El^C) illustrates the percentage of the cumulative magnitude of 
the harmonics of H(<f>) as a function of the choice of the highest harmonic component to be considered. The cutoff 
harmonic is given in terms of the feedback order to allow different orders of feedback to be compared to one another. 
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FIG. 3: Electrochemical experiments. (A) Time series of a single nearly-smooth oscillator (Vb = 1.110 V, R p = 650 fi). (B) 
Percentage of the first seven harmonics within the waveform of the rhythmic element. (C) The percentage of H(cj>) contained 
within the first k harmonics of H(cj>) as a function of the choice of k. The values of k have been re-centered by the feedback 
order n. (D-F) Experimental measurements of interaction functions corresponding to (D) 1 st order feedback (K = 0.07, ko — 



V, fci = 1, n = 0.013 rad/27r), (E) 2 nd order feedback (K = 1.6, k = -0.003 V, fci = 
(F) 3 rd order feedback (K = 18, k = 6.5 x 10~ 5 V, fci = 0, k 2 = V" 1 , fc 3 = 1 V~ 2 , 
the first 7 harmonic components within the measured interaction function using (G) 1' 
and (I) 3 rd order feedback. 



0, k 2 



1 V" 



T2 



0.013 rad/27r), and 
r 3 = 0.013 rad/27r). (G-I) Percentage of 
* order feedback, (H) 2 nd order feedback, 



Applying a 1 st order feedback signal to the experimental system produced an interaction function with a large 1 st 
order component and a relatively small 2 nd order component [Figs. OJD) and[3]JG)]. While the 1 st order harmonic 
only makes up 82% of H ((/>), the combination of the 1 st and 2 nd order harmonics account for 96% of its magnitude. 
When a 2 nd order feedback was used, it substantially reduced the 1 st order harmonic of H((f>) while increasing the 2 nd 
order harmonic [Figs. EJE) and^H)]. A small increase in the 3 rd order harmonic was observed due to anharmonicity. 
Together, these three components make up ~ 90% of the overall magnitude of H{<f>). 3 rd order feedback increases 
the ratio between the 3 rd and 1 st order harmonics of H{<f>) when compared to 1 st order feedback (0.126 for 3 rd order 
feedback vs 0.027 for 1 st order feedback). Due to strong second order harmonics in the waveform, non-trivial 2 nd and 
4 th order components were also observed in H{4>). In this case, the first four components account for ~ 98% of the 
overall magnitude of the interaction function. These results indicate that the overall shape of the of H(<j>) is largely 
composed of I th harmonics with I < n + 1, where n is the feedback order used to produce the function, in line with 
theoretical expectations. 

While the magnitude of the harmonics of H is controlled by the feedback order and their associated feedback gains 
{k n }, the ratio between the symmetric and anti-symmetric components of H is controlled by the feedback delay {t„}. 
This indicates that increasing the feedback delay is equivalent to shifting the phase of the corresponding components 
of the interaction function. To validate this claim experimentally, a series of interaction functions was measured using 
a two oscillator system with global first order feedback over a range of feedback delay t\ from 0.013 to 1 rad/27r. 

The base interaction function (ti = 0.013 rad/2-7r) is illustrated in Fig. HJA). As the feedback delay was increased, 
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FIG. 4: (Color online) (A) Experimental measurements (dots) with a Fourier fit (line) of an interaction function obtained 
using linear feedback (Vb = 1.165 V, R p — 650 Q, K = 0.07, ko = V, fci = 1, t\ — 0.013 rad/27r). (B) Interaction functions 
obtained with n = [0.013, 0.2,0.4,0.6] rad/27r respectively. (C) Phase shift of the interaction function and maximum value of 
the correlation coefficient as a function of feedback delay. 




Time < s ) <(> (rad/2jt) 4 (rad/27t) 

FIG. 5: (A) Time series of electrode potential (Vb = 1.165 V, R v = 650 fi). (B) Response function Z{cj>) and waveform (inset) 
of a single oscillator. (C) Target (solid line), optimized (dashed line), and measured (dots) interaction function with feedback 
parameters K = 0.0494, k = -0.0526 V, fci = 8.7376, k 2 = 16.3696 V" 1 , n = 0.21 rad/27r, r 2 = 0.68 rad/2vr. 

the interaction function was observed to shift [Fig. EjB)]. To determine the phase shift of the interaction functions 
when t\ > 0.013 rad/27r, a correlation function was calculated between the base interaction function and each of the 
shifted interaction functions. The correlation function was determined by finding the correlation coefficient between 
the shifted interaction functions and the base interaction function as the phase of the base function was rotated from 
to 2ir. The phase which produced the maximum value of the correlation coefficient was taken as the experimentally 
observed phase shift. FigureHjC) indicates that the phase shift of H is directly proportional to the feedback delay with 
a proportionality constant of 1. For each measurement, the maximum value of the correlation coefficients remained 
high (> 0.98), indicating that the overall shape of the interaction function was preserved. 

Knowing the relationship between the harmonics of the interaction function and the feedback parameters, it is 
possible to engineer a feedback that produces a desired interaction function, for example, H(<p) = sin(</> — a) — r sin(20) 
where a = 1.32 and r = 0.25. 

Before the feedback parameters associated with this interaction function can be calculated, the waveform [Fig. 
EJA)] and the response function of the oscillations must be determined. The response function was calculated using 
Eq. ^ from multiple measurements of interaction functions under different feedback conditions (usually 1 st , 2 nd 
and 3 rd order feedback). Since Eq. ((9]) does not have an analytical solution for the response function, a numerical 
optimization algorithm was used to calculate the Fourier coefficients of Z{4>) [(Fig. E^B)]. Once the response function 
was known, the feedback parameters {k n } and {r n } were optimized to achieve the desired interaction function also 
using Eq. §§§ [gj]. The interaction function produced by the optimized feedback parameters was experimentally 
determined to ensure that they produce the expected function. Figure [5]JC) compares the experimentally measured 
interaction function to the interaction function predicted by Eq. © . By calculating the Fourier coefficients of the 
experimental measurements, it was determined that a = 1.350 and r = 0.242, within 3% of their target values. 
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FIG. 6: Experiments: Effects of polynomial feedback on the interaction function. (A) Time series of electrode potential during 
electrodissolution of nickel wires in sulfuric acid (Vb = 1.195 V, R p = 650 Q) (B) Response function and waveform (inset) of 
a single oscillator. (C) Interaction function optimized to produce a 1 cluster state. (K = 0.4, ko = V, fci = 1, n = 0.014) 
(D) Interaction function optimized to produce a 2 cluster state. (K = 0.0425, fco = 14.97 V, fci = —3.265, k 2 = —66.087 
V -1 , n = 0.014, T2 = 0.368) (E) Interaction function optimized to produce a 3 cluster state. (K = 0.0424, fco = 20.747 V, 
fci = -4.142, fc 2 = -72.317 V , fc 3 = 251.744 V" 2 , n = 0.014 rad/27r, r 2 = 0.32, r 3 = 0.04) (F) Interaction function optimized 
to produce a 4 cluster state. (K = 0.0839, fc = 1.787 V, fci = -5.099, fc 2 = -34.136 V" 1 , fc 3 = 196.145 V~ 2 , fc 4 = 3139.686 
V" 3 , Ti = 0.014, r 2 = 0.44, r 3 = 0.02, r 4 = 0.36). The feedback delay times are given in rad/27r. 



C. Phase Clustering Experiments 



To engineer a cluster state in the experimental system, feedback parameters were selected such that the desired 
cluster state was stabilized. Four sets of experiments were conducted to obtain balanced cluster states composed of 
one to four clusters, using a population of 64 oscillators. Since a four cluster state requires the presence of fourth 
order harmonics in the response function, the operation voltage of the system was set at 1.195 V for each experiment 
causing weakly relaxational oscillations [Figs. and^B)]. The Fourier coefficients of the waveform and response 

function can be found in Table IIIII 

As seen in the numerical simulations (Sec. |Vj), there exists a large number of equally valid target interaction functions 
which can produce the desired cluster states. No specific target function was selected; The Fourier coefficients of the 
interaction function were optimized such that the desired cluster state was uniquely (or almost uniquely) stabilized. 
Given previous numerical results, n th order feedback was used to produce an n cluster state. Since linear feedback is 
sufficient to produce the one cluster state, no optimization was necessary in this case. For the higher order cluster 
states, a set of penalties were created to describe the fitness of the interaction function based on the distance between 
its Fourier coefficients and an acceptable range of coefficients. The fitness of H{4>) was calculated using the equations 

fitness = magnitude + penalties, 



magnitude = K , 

n— 1 

7 

penalties = P n , (35) 



Pn = 



\ Bn - LB "+ UB " | for B n < LB„ or B n > UB„, 
for LB„ < B n < UB n 
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Waveform Response Fn H (1 Cluster) H (2 Cluster) H (3 Cluster) H (4 Cluster) 



n 


Even 


Odd 


Even 


Odd 


Even 


Odd 


Even 


Odd 


Even 


Odd 


Even 


Odd 


1 


-0.0710 


+0.0063 


+5.6533 


+ 15.460 


-0.0799 


+0.2211 


-0.0014 


-0.1688 


-0.0619 


-0.1289 


+0.1503 


-0.0617 


2 


-0.0399 


-0.0056 


+ 12.229 


+ 14.496 


-0.1303 


+0.0819 


+0.1475 


+0.0446 


+0.1313 


-0.0157 


-0.1384 


-0.0291 


3 


-0.0212 


-0.0054 


+9.9103 


+5.7058 


-0.0507 


+0.0010 


-0.0473 


-0.0138 


-0.0657 


+0.0353 


+0.1281 


-0.0593 


4 


-0.0133 


-0.0026 


+ 1.8530 


+2.6015 


-0.0071 


+0.0029 


-0.0101 


-0.0127 


-0.0140 


-0.0148 


-0.0049 


+0.0414 


5 


-0.0064 


-0.0005 


-0.8616 


+2.1364 


-0.0004 


+0.0031 


-0.0067 


+0.0013 


-0.0009 


-0.0057 


-0.0275 


+0.0069 


6 


-0.0038 


+0.0004 


+0.5878 


+ 1.7258 


-0.0010 


+0.0011 


-0.0033 


-0.0010 


-0.0021 


+0.0024 


+0.0130 


-0.0022 



1 cluster 


Ai 


A 2 


A 3 


A 4 


2 cluster 


Ai 


A 2 


A 3 


A 4 


M = 1 


-0.422 








M = 1 


+0.172 








M = 2 


+0.058 


-0.182 






M = 2 


-0.236 


-0.032 






M = 3 


+0.196 


+0.196 


-0.010 




M = 3 


-0.014 


-0.014 


+0.048 




M = 4 


+0.108 


+0.159 


+0.108 


-0.012 


M = 4 


-0.051 


+0.134 


-0.051 


+0.051 




3 cluster 


Ai 


A 2 


A 3 


A 4 


4 cluster 


Ai 


A 2 


A 3 


A 4 


M = 1 


+0.127 








M = 1 


+0.111 








M = 2 


+0.025 


+0.076 






M = 2 


-0.299 


-0.094 






M = 3 


-0.245 


-0.245 


-0.121 




M = 3 


+0.231 


+0.231 


+0.191 




M = 4 


+0.034 


+0.043 


+0.034 


+0.059 


M = 4 


-0.268 


-0.237 


-0.268 


-0.166 



TABLE III: (Top) Fourier coefficients of the waveform, response function, and the optimized interaction functions for each of 
the four experimental objectives. (Bottom) Transversal eigenvalues for cluster states 1-4 for each of the four experiments, as 
calculated from the Fourier coefficients of the corresponding interaction function. 



Harmonic: 




2 nd 


3 rd 


4 th 


5 th 


6 th 


nth 


2 Cluster Exp. 
2 Cluster Exp. 


Lower Bounds (LB) 
Upper Bounds (UB) 


-1.0 
-0.5 


0.5 
1.0 


-2.0 
-0.4 


-0.8 
-0.3 


-0.5 
-0.1 


-0.50 
-0.05 


-0.50 
-0.05 


3 Cluster Exp. 
3 Cluster Exp. 


Lower Bounds (LB) 
Upper Bounds (UB) 


-1.0 
-0.5 


-2.0 
-0.4 


0.5 
1.0 


-0.8 
-0.3 


-0.5 
-0.1 


-0.50 
-0.05 


-0.50 
-0.05 


4 Cluster Exp. 
4 Cluster Exp. 


Lower Bounds (LB) 
Upper Bounds (UB) 


-1.0 
-0.5 


-2.0 
-0.4 


-0.8 
-0.3 


0.2 
0.5 


-0.5 
-0.1 


-0.50 
-0.05 


-0.50 
-0.05 



TABLE IV: Range of the odd Fourier coefficients of H{<j>) used to optimize feedback parameters to produce phase cluster states 
2-4. 

where Bi is the i th odd Fourier coefficient of H((f>). The upper and lower bounds (UB and LB) of the odd Fourier 
coefficients were selected such that the desired cluster state would be stable and the other (up to 6) cluster states 
in the system would be unstable. As previously demonstrated, this requires that the interaction function for an n 
cluster state have a large positive n th order harmonic, and sufficiently negative m th harmonics (m ^ n) to destabilize 
all other cluster states. The target Fourier coefficient ranges reflect this requirement, and are tabulated in Table 
IIVI Additionally, a magnitude adjustment was added to penalize parameter sets which produced a large amplitude 
feedback signal. Large feedback perturbations are not desirable since they may change the amplitude of the rhythmic 
elements of the system, violating the phase approximation. By minimizing the value of the fitness variable, the 
optimization forced the interaction function to have Fourier coefficients necessary to produce the desired cluster state. 
The optimized interaction functions arc illustrated in Figs. EIC)-[BUF). 

The transversal eigenvalues of states with one to four clusters can be seen in Table [TTTI for each experiment. They 
were calculated from the Fourier coefficients of the experimental interaction functions using Eqs. (|A3[) and (|A4|) . The 
eigenvalues indicate that the desired one, two, and three cluster states are uniquely stable. In the case of the four 
cluster experiment, the numerical optimization was unable to find feedback parameters to stabilize the four cluster 
state without also stabilizing the two cluster state. This is not unexpected, given that the difficulty of the optimization 
dramatically increases with feedback order. Therefore, the four cluster experiment will have a bi-stability between 
the four cluster state and the two cluster state. In this case, the final state of the system will be determined by the 
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FIG. 7: (Color online) (A) Time series of the Ri order parameter, using feedback optimized to produce a one cluster state. 
Arrows indicate the application and termination of the feedback signal. (B) Time series of the R2 order parameter using 
feedback optimized to produce a two cluster state. (C) Time series of the -R3 order parameter using feedback optimized to 
produce a three cluster state. (D) Time series of the R4, order parameter using feedback optimized to produce a four cluster 
state. 

initial conditions of the system. 

After the feedback parameters were determined, they were applied to the experimental system, driving it towards 
the appropriate cluster state (Fig. [?])■ Initially no feedback is present in the system, and the rhythmic elements 
were isolated from one another. Without feedback, these elements have a base frequency of 0.5 Hz ±5% with phases 
randomly distributed between and 2tt. Upon application of the feedback signal, the system progresses towards the 
desired cluster state after a short transient period. When the feedback is removed, the system relaxes back to its 
original unstructured configuration. Each experiment was successful in producing the desired cluster state from the 
appropriate feedback signal. It is important to note that although the four cluster experiment was predicted to have 
a bistability between the two and four cluster states, only the four cluster state was expcricntially observed. This 
seems to indicate that the basin of attraction for the four cluster state is sufficiently larger than the basin of the the 
two cluster state. 



VII. CONCLUDING REMARKS 

We have presented a framework for engineering target dynamical behavior in populations of oscillators with mild 
feedback. Using a time delayed, nonlinear feedback, Eq. ([5]), a variety of collective dynamics possible in weakly coupled 
oscillators can be engineered. The comprehensive theory, based on phase models, behind the methodology has been 
presented. We have verified the theoretical arguments by both numerical and experimental studies, showing that the 
methodology can be applied accurately to limit-cycle oscillator systems. As an illustration, by introducing the global 
feedback given as Eq. @, various clustering behaviors have been demonstrated numerically and experimentally 

Our methodology is based on the fact that the existence and stability conditions of dynamical states in weakly 
coupled identical oscillators are characterized by the phase interaction function. Thus, knowing an interaction function 
resulting in a target collective dynamics, the only remaining issue is how to construct the physical interaction yielding 
the phase interaction function. An interaction function can be constructed using the proposed feedback function, Eq. 
([5]) . The choice of the specific form of the feedback function was motivated by the flexible application of the imposed 
interaction function for synchronization engineering (Sec. HV[) . It has been shown that the n th order term of the 
feedback signal effectively enhances the n th Fourier components of the interaction function. The time-delay of the n th 
term is utilized to arbitrarily tune the balance of even and odd parts of the Fourier components. These correspondences 
appear intuitively reasonable, as the n power of the harmonic signal makes a component of harmonic signal having 
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n times frequency and the time delay shifts the waveform. In general, the higher order harmonics in the interaction 
function are responsible for complex dynamical behavior including dynamical clustering. Our methodology provides 
a framework for tuning all the harmonics in the interaction function. 

A major advantage of our methodology is that the feedback resulting in a target interaction function can be designed 
through the knowledge of the macroscopic observables of an isolated oscillator, that is, the waveform and the phase 
response function. When focusing on engineering synchronization properties, a microscopic investigation of the system 
is not needed. This point is beneficial when applications to biological systems are considered. It is usually a formidable 
task to construct an appropriate, detailed mathematical model of a biological system. However, the investigation of 
the phase response function is often possible; the PRC's of circadian oscillators with respect to light or temperature 
stimuli have been extensively measured [29[ as well as the PRC's of neurons with respect to electric stimuli [2 ll. [23|. 

Our methodology may be used not only to induce dynamical order but also to destroy synchronization. In a previous 
paper Q, we have demonstrated that a theoretically designed feedback successfully desynchronizes a population of 
chemical oscillators which otherwise shows simple synchronization due to global interaction among elements. The 
model-engineered feedback may find application in pacemaker and anti-pacemaker design for medical use (tremors, 
epilepsy) . 

Because of the robustness of phase description of limit-cycle oscillators, our methodology for designing interaction 
functions with feedback is robust against (at least weak) noise. However, when a complex dynamical structure is 
designed, we need to consider the (structural) stability of the designed dynamical behavior in the presence of noise. 
For example, global noise can enhance the extent of phase synchronization [30j . but can destroy subtle structures like 
slow switching [H,[l{|. Therefore, the precision of the fitted interaction function and the overall gain shall be carefully 
chosen in the presence of noise to obtain the desired structure. Because the proposed methodology was shown to work 
in the experimental system, our method should be applicable in systems with weak noise and well-defined oscillator 
waveform and response function. 

Limitations to our approach should be noted. We have focused on mild engineering, mild such that essential 
dynamical properties of elements are preserved. This strategy allows us to use the phase model. The phase models 
cannot be used with strong feedback because of amplitude effects. Also, the applicability of our method to chaotic 
oscillators is unclear because the rigorous phase description for chaotic oscillators has not been established yet. In 
coupled chaotic oscillators, various types of collective behavior arise and some of them are analogous to those in 
weakly coupled limit-cycle oscillators, such as phase synchronization [3III32I]. It would be thus worth trying to extend 
the method to chaotic oscillators. Another issue arises in cases where limit-cycle oscillators have inherent complex 
interactions. In the present paper, oscillators are assumed to be independent (i.e., uncoupled) unless feedback is 
applied. In the presence of inherent global coupling, we have shown the desynchronization is possible using our 
methodology Q . What happens if the oscillators are coupled via space dependent interactions or complex networks? 
This issue requires further exploration, for example, in chemical reaction-diffusion systems (33j and control neural 
networks. 
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APPENDIX A: EXISTENCE AND STABILITY OF THE BALANCED CLUSTER STATES 

The balanced n cluster state may be described as 



where the set Ck identifies the oscillators forming the cluster k (k = 0, . . . , n — 1) and each set includes N/n elements. 
Such a solution always exists in the phase model Substituting (|A1[) into Eq. we obtain 
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The linear stability problem for the balanced cluster states has been studied by Okuda [lTf. The eigenvalues were 
found to be 

oo 

AtL = -2E M * ( A3 ) 
1=1 



RcA intcr , p = A|" t } ra - ^2 llm { H n(i-i)+ P + H nl _ p }, (A4) 
i=i 



Ao = (A5) 

where A;™^ is associated with intra-cluster fluctuations (N — n multiplicity), A; n t or . p (p = 1, . . . , n — 1) are associated 
with inter-cluster fluctuations, and Ao is associated with the identical phase shift. 

For the interaction function with ImH n > and ImiT; < for I ^ n, the following relation holds: y 

RcA intcr , p < A^ ra = A[ n ") x . (A6) 

(n) 

Thus, the n cluster state is linearly stable if and only if Amax < 0. 



[i] 

[2] 
[3] 
[4] 
[5] 
[6] 
[7] 

[8] 
[9] 

[io; 
[ii 

[12 
[13 
[14 

[15 

[16 

[17 
[18 
[19 
[20 
[21 
[22 
[23 
[24 

[25 
[26 
[27 
[28 
[29 
[30 
[31 
[32 



A. T. Winfree, The Geometry of Biological Time (Springer, New York, 1980). 

P. A. Tass, Phase Resetting in Medicine and Biology: Stochastic Modeling and Data Analysis (Springer, Berlin, 1999). 
T. Carroll, J. Heagy, and L. Pecora, Phys. Rev. E 54, 4676 (1996). 

B. Blasius, A. Huppert, and L. Stone, Nature 399, 6734 (1999). 

Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence (Springer, New York, 1984). 

G. Ertl, Science 254, 1750 (1991). 

I. R. Epstein and J. A. Pojman, Chemical Dynamics: Oscillations, Waves, Patterns, and Chaos (Oxford Univ. Press, 
Pxford, 1998). 

A. S. Mikhailov and K. Showalter, Phys. Rep. 425, 79 (2006). 
I. Z. Kiss, C. G. Rusin, H. Kori, and J. L. Hudson, Science 316, 1886 (2007). 
Z. Neda, E. Ravasz, Y. Brechet, T. Vicsek, and A. Barabasi, Nature 403, 849 (2000). 
S. M. Reppert and D. R. Weaver, Nature 418, 935 (2002). 

S. H. Strogatz, D. M. Abrams, A. McRobie, B. Eckhardt, and E. Ott, Nature 438, 43 (2005). 
V. K. Vanag, L. Yang, M. Dolnik, A. M. Zhabotinsky, and I. R. Epstein, Nature 406, 389 (2000). 

M. Kim, M. Bertram, M. Pollman, A. von Oertzen, A. S. Mikhailov, H. H. Rotermund, and G. Ertl, Science 292, 1357 
(2001). 

S. C. Manrubia, A. S. Mikhailov, and D. H. Zanette, Emergence of Dynamical Order: Synhronization Phenomena in 
Complex Systems (World Scientific Sigapore, Sigapore, 2004). 
D. J. Christini and L. Glass, Chaos 12, 732 (2002). 
K. Okuda, Physica D 63, 424 (1993). 

D. Hansel, G. Mato, and C. Meunier, Phys. Rev. E 48, 347 (1993). 

H. Kori and Y. Kuramoto, Phys. Rev. E 63, 046214 (2001). 

H. Kori, Phys. Rev. E 68, 021919 (2003). 

R. F. Galan, G. B. Ermentrout, and N. N. Urban, Phys. Rev. Lett. 94, 158101 (2005). 

I. Z. Kiss, Y. M. Zhai, and J. L. Hudson, Phys. Rev. Lett. 94, 248301 (2005). 

Y. Tsubo, M. Takada, A. Reyes, and T. Fukai, European Journal of Neuroscience 25, 3429 (2007). 

E. M. Izhikevich, Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting (MIT Press, Cambridge, 
2007). 

E. M. Izhikevich, Phys. Rev. E 58, 905 (1998). 

P. Glansdorff and I. Prigogine, Thermodynamic theory of structure, stability, and fluctuations (Wiley, London, 1971). 
D. Golomb, D. Hansel, B. Shraiman, and H. Y. Sompolinsky, Phys. Rev. A 45, 3516 (1992). 
J. Miyazaki and S. Kinoshita, Phys. Rev. Lett. 96, 194101 (2006). 

C. H. Johnson, Chronobiology international 16, 711 (1999). 

C. Zhou, J. Kurths, I. Kiss, , and J. Hudson, Phys. Rev. Lett. 89, 014101 (2002). 
S. Boccaletti, J. Kurths, G. Osipov, D. Valladares, and C. Zhou, Phys. Rep. 366, 1 (2002). 

A. Pikovsky, M. Rosenblum, and J. Kurths, Synchronization: A universal concept of nonlinear sciences (Cambridge 
University Press, Cambridge, 2001). 
[33] Y. Kobayashi and H. Kori, submitted 



